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RESEARCH  OBJECTIVES 


The  objective  of  this  work  is  to  develop  a  novel  technique  for  studying  transient 
separated  flowfields,  such  as  those  typical  for  maneuvering  aircraft.  Transient  aircraft 
maneuvers  may  result  in  augmented  aerodynamic  forces,  yet  difficulties  in  the  prediction 
and  control  of  these  unsteady  forces  prevent  them  from  being  practically  utilized.  The 
ultimate  goal  of  this  research  is  to  develop  flow  control  techniques  using  the  boundary 
conditions  in  a  Navier-Stokes  calculation  as  a  guide. 

In  performing  a  Navier-Stokes  calculation,  the  interior  flow  is  determined  by  the 
boundary  conditions  specified  at  the  edges  of  a  domain  (neglecting  uncertainties  such  as  the 
turbulence  model,  artificial  dissipation  effects,  etc.).  As  only  a  few  parameters  can  be 
specified  on  the  boundaries,  the  variation  of  these  conditions  should  allow  the  computation 
of  any  physical  flowfield  in  the  interior.  A  method  of  determining  the  boundary  conditions 
which  are  capable  of  controlling  unsteady  separated  flow  is  sought.  Once  numerical 
boundary  conditions  are  established,  then  their  physical  counterparts  may  be  found. 

One  objective  is  to  compute  the  flow  around  a  pitching,  rolling  wing  to  use  as  a 
secondary  base  on  which  to  apply  the  flow  control  techniques  described  below.  The 
primary  base  is  a  two-dimensional  computation  which  is  less  expensive  to  use,  but  which 
contains  a  sufficient  degree  of  flow  complexity. 

There  are  two  further  objectives  in  this  work.  The  second  objective  is  the  development 
of  a  subdomain  technique  in  a  computational  investigation.  A  subdomain  of  a  flowfield  is 
chosen  such  that  it  encompasses  a  flow  phenomenon  of  interest.  The  boundary  conditions 
of  the  subdomain  allow  a  localized  study  of  the  effect  of  a  certain  flow  parameter  on  the 
control  of  this  portion  of  the  flowfield.  Also,  the  small  grid  associated  with  a  subdomain 
(compared  with  a  global  solution  which  extends  into  the  farfield)  will  lead  to  faster  running 
times,  which  is  necessary  for  time-intensive  development  work.  An  important  aspect  of 
research  on  the  use  of  subdomains  involves  the  determination  of  appropriate  boundary 
conditions  which  will  provide  a  consistent  transition  between  a  subdomain  and  a  global 
calculation. 

The  third  objective  is  the  development  of  an  intelligent  search  technique  to  find 
innovative  boundary  conditions.  The  purpose  of  this  work  is  to  find  which  boundary 
conditions  will  control  a  flowfield.  As  it  is  inefficient  to  randomly  test  new  boundary 
conditions  in  the  hope  of  discovering  a  useful  one,  a  technique  to  methodically  identify 
promising  boundary  conditions  is  required.  A  computerized  search  method,  using  artificial 
intelligence  techniques,  will  search  a  field  of  possible  boundary  conditions  in  the  most 
efficient  manner.  In  addition,  this  computerized  search  technique  may  be  standardized  to 
be  used  in  any  flow  situation. 

STATUS  OF  RESEARCH  EFFORT 

Unsteady  Flow  Computations 

The  computation  of  the  flow  around  a  pitching,  rolling  wing  was  not  achieved.  This 
was  due  to  two  problems;  first  the  construction  of  the  correct  boundary  conditions  for  the 
subdomains  proved  much  more  difficult  than  expected,  and  second,  there  was  considerable 
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uncertainty  as  to  the  availability  of  the  computer  resources  to  achieve  this  objective  to 
adequate  accuracy.  The  study  was  therefore  restricted  to  the  two-dimensional  unsteady 
flow  around  an  airfoil. 

Subdomain  Technique 

The  determination  of  appropriate  boundary  conditions  for  a  subdomain  calculation  is 
an  important  aspect  of  this  work.  The  flowfield  at  a  subdomain  boundary  is  generally  two- 
or  three-dimensional,  and  large  gradients  of  velocity  or  thermodynamic  quantities  may 
exist.  The  widely-used1  boundary  condition  theory  for  Navier-Stokes  codes  was  developed 
for  one-dimensional  inviscid  flow,  and  it  is  usually  adequate  for  far-field  boundary 
conditions  when  the  flow  is  nearly  one-dimensional.  However,  a  multi-dimensional 
boundary  condition  theory  is  needed  for  a  subdomain  calculation.  In  addition,  the 
boundary  conditions  for  a  subdomain  must  somehow  be  coupled  with  the  flowfield  within 
the  subdomain.  That  is,  if  an  upstream  or  airfoil  surface  boundary  condition  is  altered  in 
order  to  change  the  internal  flowfield,  the  downstream  boundary  condition  must  also  reflect 
the  internal  changes  so  as  not  to  be  physically  inconsistent. 

The  boundary  conditions  most  commonly  used  for  Navier-Stokes  computations  are 
derived1  from  the  one-dimensional  inviscid  gas  dynamic  equations.  These  boundary 
conditions  are: 

1.  inflow:  specify  p,  U,  V  and  extrapolate  p  from  the  interior. 

2.  outflow:  specify  p  and  extrapolate  p,  U,  and  V  from  the  interior. 

Where  U  and  V  are  the  velocity  components,  p  is  the  density  and  p  is  the  pressure. 

An  analysis  of  boundary  conditions  for  the  two-dimensional  steady  Euler  and  Navier- 
Stokes  equations  has  been  completed,2  This  note  is  attached  as  Appendix  A.  The  boundary 
conditions  for  the  steady  two-dimensional  Euler  equations  are  for  an  inflow  boundary: 

1.  pU  or  pV  is  specified 

2.  Poo  and  Poo  are  specified  (p^  p^  are  datum  values) 

3.  entropy  is  specified 

At  an  outflow  boundary,  the  entropy  is  extrapolated  from  the  interior,  rather  than  specified 
independently;  otherwise  the  outflow  boundary  conditions  are  the  same  as  those  for  inflow. 
For  one-dimensional  flow  these  boundary  conditions  reduce  to  the  conditions  given  by 
Yee.1  In  Reference  2  the  conditions  are  given  for  which  one-dimensional  boundary 
conditions  are  adequate.  All  quantities  that  are  not  specified  are  extrapolated  from  the 
interior  solution. 

The  two-dimensional  steady  Navier-Stokes  equations  may  be  written  in  a  form  to  show 
that  the  appropriate  boundary  conditions  are: 

1.  pU  and  pV  specified  on  a  boundary 

2.  either  entropy,  pressure,  or  temperature  specified  on  a  boundary  (allowing  the 
computation  of  p),  and 
both  pOT  and  p^  specified  as  a  datum. 


3. 


These  conditions  differ  from  those  given  in  Reference  1  because  the  latter  are  derived  from 
the  one-dimensional  inviscid  Euier  equations. 

Having  developed  the  boundary  conditions  it  becomes  a  more  than  trivial  matter  to 
numerically  implement  them.  In  particular  the  correct  boundary  conditions  for  a 
subdomain  are  difficult  to  implement  because  it  is  almost  impossible  to  specify  values  of 
pU,  pV,  and  p  on  the  subdomain  boundary  that  are  consistent  with  the  flow  in  the  adjoining 
subdomain.  A  wide  range  of  existing  boundary  condition  treatments  were  tested  with  little 
success;  these  techniques  and  extensions  of  these  techniques  are  noted  below. 

Computations  over  a  subdomain  of  a  global  solution  field  have  shown  that  the  interior 
flow  solution  for  the  subdomain  is  critically  dependent  on  the  outflow  boundary  conditions. 
For  example,  if  the  flow  over  the  airfoil  is  separated,  but  the  outflow  boundary  condition  is 
attached  flow,  the  interior  flow  will  reattach  before  reaching  the  back  boundary.  This 
contrasts  with  the  global  solution,  where  the  flow  remains  separated.  A  subdomain 
calculation  such  as  the  one  just  described  will  probably  be  unstable  and  will  diverge  after 
some  time.  All  of  the  calculations  that  were  performed  with  outflow  boundary  conditions 
that  were  physically  inconsistent  with  the  interior  flow  were  unstable. 

It  is  possible  to  reproduce  a  global  solution  in  a  subdomain  if  data  from  the  global 
solution  is  used  as  the  boundary  conditions.  However,  using  all  of  the  global  data  that  is 
available  will  overspecify  the  boundary  conditions.  If  partial  data  is  used  for  the  boundary 
conditions  instead,  such  as  specifying  three  variables  at  an  inflow  and  one  at  an  outflow,  the 
interior  flow  will  be  mostly  reproduced,  with  the  exception  of  the  outflow  regions.  Those 
regions  will  be  somewhat  different  from  the  global  solution.  If  the  boundary  is  in  a  region 
of  high  flow  gradients,  the  more  freedom  that  is  allowed  in  specifying  the  boundary 
conditions,  the  more  the  boundary  conditions  will  differ  from  the  global  data  they  were 
derived  from.  For  instance,  if  the  global  data  is  used  to  provide  boundary  conditions  in  the 
form  of  the  Riemann  invariants,  the  resultant  flow  and  boundary  conditions  will  have 
enough  flexibility  to  change  significantly  from  the  original  global  solution. 

Boundary  conditions  applied  in  a  wake  region  were  investigated  for  use  on  a 
subdomain  boundary.  It  was  thought  that  formulating  the  boundary  conditions  in  terms  of 
the  Riemann  invariants  would  produce  a  reasonable  result,  since  there  is  more  flexibility  in 
that  formulation  to  allow  the  flow  to  adjust  to  the  interior.  One  problem  with  this  technique 
is  that  the  Riemann  invariants  are  derived  for  one-dimensional  flow,  and  there  is  no 
corresponding  formulation  for  two-dimensional  flow.  Consequently,  there  is  some  queshon 
as  to  what  pseudo-one-dimensional  velocity  to  choose  in  the  formulation  of  the  invariants. 
Either  the  normal  velocity  to  the  boundary  or  the  total  velocity  seems  to  be  the  best  choice, 
although  there  is  no  strict  justification  for  either. 


The  formulation  of  these  boundary  conditions  that  is  used  in  the  NASA  Ames  Navier- 
Stokes  code  ARC2D  is 


Rl  =  Qn  -  2a/(y  - 1) 
R2  =  Qn  +  2a/(y  - 1 ) 
entropy 

Q, 


specified  at  inflow 
specified  at  outflow 
specified  at  inflow 
specified  at  inflow 


where  Qn  and  Qt  are  the  normal  and  tangential  velocity  components  relative  to  the 
boundary. 

This  formulation  was  tried  using  either  Qn  or  Qt  t  ,  in  the  equations  for  the  Riemann 
invariants.  It  was  found  that  in  a  wake,  where  there  are  large  flow  gradients  present  in  the 
vicinity  of  the  boundary,  the  changes  in  flow  quantities  between  the  first  and  second  grid 
point  were  large  enough  to  cause  very  large  errors  in  the  boundary  conditions.  For 
instance,  for  an  inflow,  Rl  is  specified  at  the  boundary  and  R2  is  extrapolated  from  one  grid 
point  away  from  the  boundary.  Rl  and  R2  are  combined  to  solve  for  the  speed  of  sound,  a. 
Since  the  speed  of  sound  is  close  to  one,  and  the  density  goes  as  (a2  )l  1  (  V' 1  * ,  a  2%  error  in 
the  speed  of  sound  leads  to  a  10%  error  in  the  density.  Adding  constraints,  such  as  the 
removal  of  the  dependence  of  the  speed  of  sound  on  the  Riemann  invariants,  improved  the 
solution  somewhat  but  still  did  not  give  satisfactory  results. 

From  the  analysis  in  Reference  2,  certain  flow  quantities  or  their  derivatives  may  be 
specified  as  boundary  conditions.  The  boundary  conditions  can  consist,  therefore,  of  either 
setting  a  flow  quantity,  its  first  derivative,  or  its  second  derivative  equal  to  zero  on  the 
boundary.  These  boundary  conditions  represent  a  case  where  all  information  is  taken  from 
the  interior  flow,  regardless  of  whether  there  is  inflow  or  outflow  across  the  boundary. 
These  boundary  conditions  were  generally  stable,  but  did  not  result  in  solutions  which 
resembled  the  global  solution. 

An  attempt  was  made  to  formulate  the  boundary  conditions  along  a  wake  in  terms  of 
the  vorticity  and  entropy,  as  these  quantities  are  conserved  along  streamlines  in  two- 
dimensional  inviscid  flow.  It  was  felt  that  using  vorticity  as  a  boundary  condition  would  be 
convenient  for  flow  control  problems,  since  vortex  generators  are  commonly  used  for  flow 
control.  The  Riemann  invariants  were  also  needed  to  complete  the  boundary  condition  set. 
Unfortunately,  since  the  vorticity  must  be  integrated  to  find  the  velocities  and  the  Riemann 
invariants  only  provide  information  on  the  total  velocity,  information  about  the  sign  of  one 
of  the  velocity  components  is  lost.  Another  problem  with  this  formulation  is  that  it 
appeared  to  be  unstable  for  any  slight  perturbations  in  the  solution,  and  it  was  thus 
unusable  for  a  subdomain  boundary  condition. 

Nonreflecting  boundary  conditions  were  investigated  for  use  on  the  outflow  boundary. 
These  boundary  conditions  suppress  reflections  from  outgoing  waves  and  cancel  incoming 
waves  for  a  domain.  Nonreflecting  boundary  conditions  have  been  derived  in  the  literature 
for  hyperbolic  equations  on  cartesian  grids.3  In  this  work,  the  boundary  conditions  were 
derived  for  use  on  curvilinear  grids  and  were  applied  to  low  speed  vortex  shedding 
problems.4'5  Appendices  B  and  C  describe  the  results  from  this  work.  Some  of  these 
results  were  also  presented  at  the  Minisymposium  on  Outflow  Boundary  Conditions, 
Swansea,  United  Kingdom,  July  10, 1989. 

The  trials  with  nonreflecting  boundary  conditions  showed  that  these  conditions  are 
inappropriate  for  subsonic  vortex  shedding  flows.  For  a  single  vortex  passing  through  a 
boundary,  the  convected  quantities  are  computed  satisfactorily,  but  the  pressure  is  highly 
inaccurate.  For  a  vortex  street,  the  frequency  of  the  vortex  shedding  is  rapidly  changed 
from  the  global  solution. 


The  primary  reason  for  the  failure  of  the  boundary  conditions  described  above  is  that 
they  fail  to  estimate  adequately  the  flow  outside  the  computational  subdomain. 
Nonreflecting  boundary  conditions  simply  allow  the  perfect  transmission  of  waves  through 
the  boundary  but  do  not  consider  whether  the  boundary  condition  implies  a  physically 
unreasonable  solution  outside  the  computational  domain.  For  certain  conditions,  a 
boundary  condition  that  removes  this  uncertainty  was  derived  during  the  present  work. 
The  idea  is  described  briefly  below;  a  more  complete  description  is  given  in  Appendix  D. 

If  a  particular  flow  is  started  from  zero,  then  the  flow  in  the  exterior  domain  is  exact 
until  waves  from  the  disturbance  in  the  interior  domain  cross  the  boundary.  If  true 
nonreflecting  boundary  conditions  could  be  devised,  then  the  outflow  boundary  condition 
would  be  correct  if  the  computation  was  started  at  zero  time  and  is  time  accurate.  An  error 
in  the  boundary  conditions  due  to  approximations  in  the  formulation  can  give  rise  to  a 
nonphysical  flow  in  the  exterior  domain  and  hence,  ultimately,  the  interior  domain.  A 
principal  difficulty  is  in  finding  an  analytic  relationship  that  will  give  the  boundary 
conditions  for  the  Navier-Stokes  equations.  The  progress  in  nonreflecting  boundary 
conditions  shows  that  finding  a  true  two-dimensional  boundary  condition  is  beyond  the 
present  state  of  the  art. 

The  new  idea  is  based  on  the  following  facts: 

(a)  If  the  solution  is  started  from  a  zero  disturbance  state,  then  the  outflow  boundary 
condition  is  correct  until  the  disturbance  waves  pass  the  boundary. 

(b)  The  flow  in  the  exterior  domain  is  determined  solely  by  its  inflow  boundary 
conditions  which  are  compatible  with  the  outflow  boundary  condition  for  the 
interior  domain. 

It  is  also  assumed  that  the  flow  in  the  exterior  domain  can  be  constructed  by  a 
superposition  of  flow  with  linearized  inflow  boundary  conditions.  The  range  of  validity  of 
this  assumption  is  not  clear  but  the  present  results  are  very  encouraging. 

Based  on  these  assumptions,  the  flow  in  the  exterior  domain  is  determined  at  each  time 
step  by  the  flow  parameters  on  the  penultimate  grid  line  in  the  computational  domain.  The 
only  part  of  the  exterior  flow  that  is  important  is  that  on  the  outflow  boundary  of  the 
interior  solution  which  is  assumed  to  overlap  with  the  exterior  domain.  The  exterior  flow  is 
computed  by  relating  the  time  dependent  flow  to  a  superposition  of  the  flow  due  to  a  step 
change  on  the  boundary  (the  indicial  responses).  Results  for  the  unsteady  transonic  small 
disturbance  equation  for  a  pitching  airfoil  are  given  in  Appendix  D. 

Flow  Control  and  Optimization 

Attempts  were  made  to  steady  the  vortex  shedding  process  behind  a  stalled  airfoil,  in 
order  to  hold  a  vortex  fixed  over  the  upper  surface.  If  a  vortex  could  be  held  fixed  or 
slowed  down,  the  high  lift  over  an  airfoil  could  be  controlled.  Various  combinations  of  jets 
or  suction  over  the  upper  surface,  both  normal  and  tangential,  were  tried.  It  was  also 
attempted  to  hold  the  vortex  steady  by  fixing  the  instantaneous  boundary  conditions  about 
a  subdomain  containing  the  vortex.  None  of  these  attempts  worked.  It  does  not  seem  to  be 
possible  to  use  steady  flow-through  boundary  conditions  for  an  unsteady  flow,  even  if  it  is 


desired  to  have  the  resultant  solution  be  steady.  Especially,  it  is  not  possible  to  use  an 
instantaneous  solution  to  an  unsteady  flow  and  for  the  boundary  conditions  to  a  steady 
flow.  Since  a  solution  to  the  unsteady  Navier-Stokes  equations  is  not  necessarily  a  solution 
to  the  steady  Navier-Stokes  equations,  then  if  a  steady  vortex  flow  can  be  found,  it  will 
probably  be  a  somewhat  different  flowfield  than  what  was  started  with.  Flexible  boundary 
conditions  are  necessary  for  consistency  with  the  interior  solution,  and  this  leads  to  the 
problem  of  boundary  c.  ndition  specification  again. 

A  study  of  optimization  techniques  was  made  to  search  for  boundary  conditions  which 
have  a  controlling  effect  on  a  flowfield.7  This  study  was  made  independently  of  the 
outflow  boundary  condition  work.  Since  efficiency  is  a  major  concern  when  faced  with  the 
possibility  of  many  iterations  of  a  CFD  solution,  it  was  decided  to  use  artificial  intelligence 
search  techniques  rather  than  a  purely  mathematical  optimizer. 

A  basic  study  with  a  mathematical  optimizer  was  undertaken  initially,  in  order  to  show 
that  optimizing  techniques  may  be  used  to  find  desirable  sets  of  boundary  conditions.  The 
test  problem  was  that  of  blowing  over  the  upper  rear  surface  of  an  airfoil.  The  magnitude  of 
the  blowing  and  the  area  were  allowed  to  vary,  and  the  objective  was  to  maximize  the  lift 
coefficient  and  minimize  the  jet  area.  Although  this  was  a  simple  test  case,  the  results 
confirmed  that  mathematical  optimizers  are  inefficient  and  tend  to  get  stuck  on  local  optima 
rather  than  on  the  global  optimum. 

The  decision  to  use  an  artificial  intelligence  search  technique  was  made  because  it  could 
be  programmed  to  give  information  about  the  various  dependencies  of  the  solution  on  the 
input  variables,  as  well  as  the  optimal  values  of  the  inputs.  The  technique  that  was 
developed  here  acts  as  a  postprocessor  on  a  series  of  CFD  results.  Each  computation  uses  a 
slightly  different  input,  and  various  results  for  each  input  are  documented.  The 
postprocessing  technique  then  finds  the  dependencies  of  the  results  on  each  input  and  how 
the  combination  of  inputs  changes  the  solution. 

This  technique  was  tested  on  the  flow  control  of  an  airfoil  with  unsteady  separation. 
The  airfoil  has  pulsed  blowing  or  suction  over  the  upper  rear  of  the  airfoil.  The  inputs  that 
were  varied  with  each  calculation  were  the  magnitude  of  the  blowing  or  suction,  the 
blowing  pulse  duration,  and  the  interval  between  pulses.  The  objective  was  to  increase  the 
average  forces  on  the  airfoil  and  to  decrease  the  rms  forces.  This  technique  found  an 
optimum  for  a  low  suction  magnitude  and  short  duration  and  interval  pulses.  It  also 
identified  that  the  solution  was  most  sensitive  to  the  inputs  for  long  blowing  durations, 
short  intervals  between  blowing  pulses,  and  high  suction  values. 

A  detailed  description  of  this  method  is  found  in  Appendix  E.  This  technique  may  be 
applied  to  any  complex  flow  field  and  is  not  limited  by  the  number  of  variables  in  the 
problem.  The  possibility  exists  of  extending  this  technique  to  finding  the  interrelationships 
among  several  time-varying  flow  quantities,  rather  than  between  the  output  quantities  and 
the  input.  This  extension  to  the  technique  would  find  useful  applications  in  areas  such  as 
turbulence  modeling. 
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Summary 


The  main  goal  of  the  research  reported  here  is  to  develop  a  method  of  flow  control  for 
highly  maneuverable  aircraft.  The  method  used  in  the  work  is  based  on  the  idea  that  all 
flow  control  devices  are  represented  by  boundary  conditions,  and,  therefore,  only  the  small 
set  of  boundary  conditions  need  be  manipulated  to  achieve  control  rather  than  a  large  set  of 
possible  physical  control  devices.  In  order  to  determine  the  effect  of  a  control  on  an  element 
of  the  flow,  for  example,  a  vortex,  the  flow  domain  was  divided  into  subdomains,  each  of 
which  contained  one  or  more  flow  features  of  interest.  An  artificial  intelligence  optimizer 
was  developed  to  determine  the  best  controls  to  use. 

The  work  performed  during  the  contract  did  achieve  most  of  the  tasks,  such  as  the 
development  of  an  intelligence  optimizer  and  some  insight  into  methods  of  control.  The 
main  stumbling  block  was  the  inadequate  status  of  boundary  conditions  for  the 
subdomains,  a  fact  that  was  foreseen  at  the  outset.  None  of  the  currently  used  boundary 
conditions  in  CFD  proved  to  be  satisfactory  and  only  towards  the  end  of  the  contract  did  a 
possible  solution  (the  indicial  boundary  condition)  appear.  As  an  additional  point,  the  use 
of  this  boundary  condition  could  reduce  computational  time  of  CFD  calculations 
considerably  because  of  its  capability  of  reducing  the  grid  size  substantially. 
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APPENDIX  A 


Derivation  of  the  Boundary  Conditions  for  the 
Two-Dimensional  Euler  and  Navier-Stokes  Equations 

BOUNDARY  CONDITIONS  FOR  THE  EULER  EQUATIONS 

Introduction 

This  paper  is  concerned  with  deriving  the  necessary  and  sufficient  boundary 
conditions  for  both  the  Euler  and  Navier-Stokes  equations.  In  Part  I,  the  boundary 
conditions  for  the  Euler  equations  are  developed  for  two-dimensional  flow  and  are 
compared  with  the  restrictions  implied  by  the  usual1  outflow  boundary  conditions  that  are 
based  on  one-dimensional  theory.  It  is  shown  that  the  commonly  used  characteristic 
boundary  conditions  are  not  generally  relevant  to  multi-dimensional  problems. 

In  Part  II  of  this  paper,  the  boundary  conditions  for  the  compressible,  two-dimensional 
Navier-Stokes  equations  are  derived. 

Basic  Equations 
The  vorticity  equation 


VVX=" 

can  be  written  in  terms  of  a  stream  function  f  such  that 

V'xx  +Vyy  =PW+PyU-pxV  (1) 

where  oj  is  the  vorticity  and 

Vy  =  PU  -  pJJ^  fx  =  -pV  (2) 

By  manipulation  of  the  energy  and  gas  equations,  and  using  Gibbs  relation,  the  density 
is  given  in  terms  of  the  entropy,  S,  and  the  velocity  components  U,  V  by 
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he  vorticity  is  given  by  Crocco's  theorem  as 


where  R  is  the  gas  constant 
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and 


is  the  derivative  in  the  direction  normal  to  the  streamlines.  Note  that 

a  _u  a  _v  a 

a  n  q  du  q  dx 

Along  a  streamline  the  entropy  is  constant  and  hence 
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The  pressure  is  given  by 
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Using  Equations  (3),  (4),  (5),  and  (6),  Equation  (1)  can  be  written  as 


(7) 


W  +  V 
xx  yy 


v/y-l 

p  Q  Y  exp (-S/R) 

oo 

Rq2 


as  as 

(U  -r—  —  V  — — ) 
dy  dx 


(8) 


+  p  1exp(-S/R)]  U  -  p  4-  [QV  1exp(-S/R)]V 

OO  0  y  CO  d  X 


where  Q  is  a  function  of  U  and  V  and  is  given  by 


0  .  u  ♦  <lz±L  M2  a  -  2- 

Z  oo  z 

u 


^2>1 


(9) 


CO  oo 


Equation  (8)  can  be  written  as 

V'xs  +Vyy  =Kx,y) 

where  Kx,  y)  is  the  right  hand  side  of  Equation  (8). 


(10) 


Consider  the  inviscid  transonic  problem  illustrated  in  Figure  1.  The  solution  is 
described  by  the  Euler  equations.  We  assume  that  the  airfoil  can  be  represented  by 
specifying  boundary  conditions  on  the  chord  line,  that  is,  thin  airfoil  boundary  conditions 
are  applied.  An  interesting  example  of  a  subdomain  is  as  follows.  The  subdomain  is 
defined  as  the  region  bounded  by  the  lines  AH,  BC,  CD,  and  DA. 
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Equation  (10)  can  be  written  as  an  integral  equation  using  Green's  Theorem.  Thus,  in 
the  domain  D 
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where  the  kernel  function  K  is  given  by 

K  =  ^ln[(x-£)2  +(y-r?)2]1/2 
Differentiation  of  Equation  (11)  with  respect  to  y  gives 
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If  the  mass  fluxes  are  known  on  the  boundaries  then  the  line  integrals  are  known.  These  are 
the  necessary  boundary  conditions  for  the  subproblem. 


Differentiation  of  Equation  (11)  with  respect  to  x  gives 
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From  Equations  (3)  and  (8)  it  can  be  seen  that  p  and  f  are  functions  only  of  U  and  V 
and  the  entropy  S.  The  entropy  is  given  by  Equation  (6)  and  the  velocities  just  ahead  of  the 
shock  wave.  Thus  S  is  a  function  only  of  U  and  V.  Hence,  if  pV  is  defined  on  y  =  ±  0  as  a 
boundary  condition.  Equation  (13)  and  (15)  are  equations  for  U  and  V  throughout  the  flow 
field.  Because  of  the  complexity  of  the  relations  for  p,  S,  and  f,  it  is  not  easy  to  determine  the 
number  of  possible  solutions.  However,  it  is  safe  to  say  that  if  pV  or  pU  is  specified  on  the 
boundary,  then  pU  and  pV  are  determined  completely  in  the  flow  field.  pU  is  also 
determined  on  the  boundary.  In  other  words,  U  and  V  cannot  both  be  specified  on  the 
boundary. 

If  a  constant  pressure  boundary  condition  is  used  then  Crocco's  equation  gives 
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where  R  is  the  gas  constant,  Cv  the  specific  heat  at  constant  volume,  and  AS  is  the  entropy 
due  to  the  shock  wave.  If  a  small  disturbance  approximation  is  made  such  that  only  linear 
terms  involving  v  are  retained,  then  a  study  of  the  density-velocity  relation  gives 
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and  hence  on  the  boundary,  if  p  is  constant 

f(q,rj)  =  pca  +  p^U  =  (pU)^  (18) 

On  the  boundary.  Equation  (13)  therefore  reduces  to  an  equation  for  pU  in  terms  of  pV. 


It  is  of  interest  to  assume  one-dimensional  flow  for  comparison.  In  this  case, 
conditions  are  uniform  across  DA  and  the  vorticity  is  zero.  Therefore,  f  is  zero  and  Equation 
(13)  can  have  the  solution 

pU|,  =  p (p, AS) U (p, AS) | ,  =  constant.  (19) 

From  Equation  (7),  the  one-dimensional  equation  gives 
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and  Equation  (19)  becomes  an  equation  for  the  entropy  S  in  terms  of  the  given  p.  If  S  is 
known,  U  is  found  from  Equation  (20)  and  p  is  found  from  Equation  (3).  Thus  all  of  the 
unknowns  at  the  outflow  boundary  are  found  in  terms  of  the  specified  pressure. 

Inflow  Boundary  Conditions 

For  the  above  analysis  of  boundary  conditions  to  be  applicable  the  properties  of  the 
fluid  must  be  specified.  The  most  obvious  place  to  specify  this  is  at  the  inflow  boundary. 
Thus  at  the  i  iflow  boundary  the  following  must  be  specified 


(a) 

pU  or  pV 

(PooUoo  in  the  above  context) 

(b) 

P 

(Pec  in  the  above  context) 

(c) 

P 

(p<x,  in  the  above  context) 

Outflow  Boundary  Conditions 

At  the  downstream  boundary  pV  (say)  is  specified,  and  if  p  is  a  constant,  then 
Equation  (13)  is  an  equation  for  pU  on  the  boundary.  The  entropy  is  given  by  the  internal 
solution  which  then  gives  p.  Thus  at  the  outflow  boundary  the  following  must  be  specified. 

(a)  p 

(b)  pV  (or  pU) 

The  other  quantities  are  found  from  the  solution  (i.e.,  extrapolation) 
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and  hence,  using  Equation  (3) 
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Use  of  One-Dimensional  Boundary  Conditions 

Equation  (13)  reduces  to  the  one-dimensional  equation  if  f  and  V  are  zero.  Hence  one¬ 
dimensional  boundary  conditions  may  be  applicable  if 

I  pV  t  «  1  (22) 


and  if 


I  f  I  =  I  pw  +  pn U I  «  1  (23) 

Thus  the  "crossflow"  on  the  boundary  must  be  small  and  the  vorticity  must  be  small. 
Extropolation  in  the  x  direction  cannot  be  used  if  is  not  small,  and  from  Equation  (18) 
this  implies  that 


dpo 

dx 


and 


a(p,u) 

dx 


must  be  small  if  extrapolation  is  to  be  used. 


BOUNDARY  CONDITIONS  FOR  THE  NAVIER-STOKES  EQUATIONS 

The  2-D  Navier-Stokes  equations  are 
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g  is  the  coefficient  of  viscosity  and  k  is  the  coefficient  of  thermal  conductivity. 


If  the  thermal  conductivity  is  neglected  (i.e.,  k  =  0),  then  Equation  (3)  can  be  written  as 
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Equation  (2)  can  be  written  as 
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Equation  (8)  can  be  written  as 
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Equation  (9)  can  be  written  in  terms  of  the  vorticity,  w,  as 
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Using  Gibbs  relation. 
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Equation  (10)  can  be  written  as 
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In  general. 
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since  by  definition 
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then,  on  using  the  energy  equation,  Equation  (6)  becomes 
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Note  that  the  highest  derivative  of  LL  in  Equation  (16)  is  a  first  derivative. 
The  basic  equation  for  the  stream  function  is 
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and  from  Equations  (13)  and  (16)  it  is  clear  that  the  highest  derivative  of  the  velocity  in 
Equation  (1 7)  comes  from  the  term  connected  with  pu,  through 


i  ,  Hii 

v  1  ax, 

D 


,  dr.  . 

u  '  ax .  ' 

3 


Now  from  Equations  (4)  and  (5) 
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and,  thus 
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and  hence  the  highest  derivative  terms  in  Equation  (17)  are  given  by 


U  U  V  V 

U  I  _ZZ  +  _2S£ _ _ ZZ  ) 

H  v  v  u  u  ; 


If  a  stream  function  ip  is  defined  such  that 


pU  =  V'v,  pV  =  -yx< 

then  Equation  (17)  can  be  written  as 

UV2(^y)  +  W2tyx)  =  F(yx,  ry/  p,  S,  T) 

where  F  contains  lower  derivatives  than  the  third.  Equation  (20)  can  be  written  as 


V-i ii  =  f 


where 


v  =  -vl  +  vl 


A  -  h 


Equation  (21)  is  a  Poisson  equation  and  requires  either  ff  or  the  partial  derivative  of  ^  be 
specified,  i.e., 

yx  or  is  to  be  specified 

If  W  is  specified  then 

-wl  +  w] 

is  specified  which  is  a  combination  of  pU  and  pV.  To  solve  for  \ft  the  differential  Equation 
(22)  must  be  solved;  this  requires  that  f  or  one  of  its  derivatives  be  specified.  Thus  to  solve 
Equation  (21)  two  of  the  boundary  conditions  are  that  pU  and  pV  be  specified. 

Now  from  Equations  (1)  and  (2) 


au.  ,  8r. 

pu.u.  -r-i  =  -U.  -P-  +  u.  — — — ^ 

1  1  3x  ,  10X.  1  dx, 

D  1  D 


Thus  from  Equations  (6),  (23),  and  (14) 
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Hence,  on  using  continuity  (Equation  (1)),  Equation  (25)  gives 


| —  (pU.H)  =  | —  (t.  .U.) 
dx,  r  1  dx.  1]  1 

j  J 


An  alternative  version  is 


dp u  . 

(pU.H)  =  (t.  ,u.)  +  H  -r — - 

ax .  p  j  ax.  ij  i  o  ax . 

3  1  3 


pUXH  -  Hq)  =  Tj-Uj 


so  that  if  HQ  is  the  stagnation  value  of  H,  Equation  (30)  is  consistent. 

Multiplication  of  Equation  (30)  by  U-  and  using  Equations  (4)  and  (5)  gives 
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Equation  (32)  can  be  manipulated  to  give 


2  „  au, 

12  u  aq  2  k  _ 

C  pT  +  -  pq  -  C  pT  =  Hr  U  .  -  q  - -  =  Q 
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From  Gibbs  relation  and  the  gas  law 


(34) 


Y 

(£-)  /  (^-)  =  exp  (S/Cv) 

OO  QO 

then 

y 

C  pT  =  C  £  =  — ^r-  p  (^-)  exp  (S/C  )  (35) 

P  p  R  y-1  oo  p  V 

and  Equation  (33)  gives 

y 

exp  (S/C  )  (~r)  p  (£-)  +  p  (£-)  (^  q2  -  C  T  )  =  Q  (36) 

v  v-1  oo  p  oo  p  2  p  o 

oo  oo 

which  is  an  equation  for  p  in  terms  of  Q,  S,  p^,  and  p^.  If  Uj  (and  hence  Q)  is  specified  on 
the  boundary  and  if  S  is  known,  then  p,  p,  and  T  (and  consequently  h)  can  be  found  from 
Equations  (34),  (36),  and  the  gas  equation.  In  general,  S,  p,  or  T  must  be  specified  on  the 
boundary  and  the  other  thermodynamic  variables  are  given  in  terms  of  this  variable  and  p. 
The  boundary  conditions  for  the  steady  Navier-Stokes  equations  are  therefore 

(a)  pU,  pV  specified  on  a  boundary 

(b)  S,  p,  or  T  specified  on  a  boundary 

(c)  p ^  p^  specified  as  a  datum 
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APPENDIX  B 


The  Application  of  Nonreflecting  Boundary  Conditions 
to  2-D  Unsteady  Computations  on  Curvilinear  Grids 


SUMMARY 

The  ability  to  apply  flow-through  boundary  conditions  in  places  where  the  flowfield  is 
neither  steady  nor  uniform  would  allow  CFD  computations  over  smaller  domains,  thus 
requiring  fewer  grid  points.  Recently,  nonreflecting  boundary  conditions  have  received  a 
lot  of  attention  as  a  promising  method  for  specifying  outflow  boundary  conditions  on  a 
shortened  domain.1 ' 3  Nonreflecting  boundary  conditions  have  been  derived  by 
Thompson1  for  the  Euler  equations  in  Cartesian  and  spherical  coordinates.  This  work 
extends  Thompson's  formulation  to  curvilinear  coordinates,  so  that  it  may  be  applied  to 
practical  two-dimensional  CFD  problems.  The  theory  is  valid  for  hyperbolic,  quasi  one¬ 
dimensional  flows,  but  loses  accuracy  and  stability  for  computations  of  low  subsonic  flows 
with  strong  two-dimensionality.  This  loss  of  accuracy  is  demonstrated  for  the  problem  of 
subsonic  vortex  shedding  through  a  boundary. 

INTRODUCTION 

Most  computational  fluid  dynamics  (CFD)  problems  are  solved  assuming  that  the 
solution  domain  boundaries  are  either  in  a  uniform  stream  or  at  a  solid  surface.  For 
boundaries  in  a  uniform  freestream,  which  is  at  a  known  velocity  and  angle  of  attack, 
boundary  conditions  derived  from  one-dimensional  characteristic  theory  can  be  applied  in  a 
relatively  straightforward  manner.  However,  it  is  of  practical  importance  to  be  able  to 
apply  boundary  conditions  in  locations  where  there  are  flow  gradients  present,  or  where  the 
flow  is  unsteady.  For  example,  many  flows  have  a  wake  structure  which  persists  for  a  long 
distance  downstream,  and  it  is  too  computationally  expensive  to  have  a  solution  domain 
extend  to  where  uniform  flow  conditions  may  be  applied. 

Nonreflecting  boundary  conditions  have  been  developed  as  a  promising  means  of 
applying  boundary  conditions  at  an  outflow.1 ' 3  These  boundary  conditions  prevent  the 
reflections  of  outgoing  waves  back  into  the  computational  domain  and  also  eliminate 
incoming  waves,  which  are  unknown  since  they  originate  from  outside  the  domain. 
However,  these  boundary  conditions  are  only  accurate  for  a  limited  subclass  of  flow  fields. 
For  instance,  Thompson1  noted  that  these  boundary  conditions  will  be  inaccurate  for 
problems  where  the  flow  should  be  influenced  by  incoming  waves,  and  Watson  and  Myers  2 
found  that  the  formulation  used  for  these  boundary  conditions  is  unstable  for  nonplanar 
acoustic  sources  at  an  inflow.  Other  errors  may  occur  in  viscous  subsonic  flow 
computations,  as  the  boundary  conditions  are  derived  from  the  hyperbolic  Euler  equations. 

Although  some  errors  may  be  expected,  nonreflecting  boundary  conditions  provide 
one  of  the  few  options  currently  available  for  treating  the  outflow  of  unsteady  nonuniform 
flows.  Therefore,  it  is  necessary  to  investigate  the  suitability  of  these  boundary  conditions 
for  practical  flow  problems,  to  determine  whether  any  errors  involved  are  within  a  tolerable 
range.  This  paper  investigates  the  application  of  nonreflecting  boundary  conditions  to 
subsonic  flows  with  a  wake  structure  that  passes  through  the  outflow  boundary. 


DERIVATION  OF  CURVILINEAR  BOUNDARY  CONDITIONS 


The  vast  majority  of  Euler  and  Navier-Stokes  finite  difference  computations  for 
aeronautical  applications  take  place  on  curvilinear  grids,  in  order  to  allow  for  realistic 
geometries.  The  first  stage  of  this  work  was  to  extend  Thompson's1  formulation  of 
nonreflecting  boundary  conditions  to  the  general  curvilinear  case.  This  work  was 
completed  for  boundaries  in  both  the  curvilinear  coordinate  directions. 

The  first  step  was  to  write  the  Euler  equations  in  primitive  variable  form  and  then  to 
group  together  the  spatial  derivatives  which  are  normal  to  the  computational  boundary. 
These  derivatives  will  eventually  be  rewritten  to  represent  the  nonreflecting  condition, 
while  the  spatial  derivatives  tangential  to  the  boundary  will  be  left  unaltered. 

Let  x  and  y  represent  the  Cartesian  coordinate  variables  and  £  and  q  represent  the 
curvilinear  coordinate  variables.  For  a  problem  where  the  outflow  boundary  runs  along  an 
rj=constant  gridline,  the  terms  in  the  Euler  equations  involving  rj-derivatives  are  to  be 
rewritten.  The  Euler  equations  are  placed  in  characteristic  form,  so  that  the  eigenvalues  are 
given  explicitly.  The  characteristic  terms  for  the  rj-derivatives  are  given  by: 
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let  a  = 

u  and  v  are  the  velocities  along  x  and  y,  U  and  V  are  the  contravariant  velocities  along  £ 
and  q,  and  c  is  the  speed  of  sound.  This  formulation  of  the  characteristic  equations  assumed 


quasi  one-dimensional  flow/  that  is,  the  derivatives  with  respect  to  £  and  fj  are  separable, 
and  the  cross-terms  are  neglected. 

For  a  subsonic  outflow  boundary.  Ml  represents  incoming  waves,  and  M2  -  M4 
represent  outgoing  waves.  In  the  Euler  equations,  the  rj-derivative  terms  are  written  in 
terms  of  Ml  -  M4.  Then,  at  the  boundary.  Ml  is  set  equal  to  zero,  and  M2  -  M4  are 
computed  numerically.  The  equations  are  written  in  conservative  form  and  are  solved 
together  with  the  interior  equations.  The  final  form  of  the  Euler  equations  at  the  boundary 
is 
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and  T  is  a  4  x  4  matrix  made  up  of  the  following  terms: 
Tn  =V 

T12  =  (V  +  ca)pi]x/ (2c) 

Tl3  =  (V  +  ca)pny/(2c) 

Tu  =  -(V  -  car) /(2c2) 

T21  =uV 

T22  =  pV  +  (V  +  ca)putjx/(2c) 

T23  =  (V  +  ca)put]y/(2c) 

T24  =  "(V  '  ca)u/(2c2)  +  nx 

T31 

T32  =  (V  +  ca)pvtjx/(2c)  -  (V  -  ca)pt]x/(2rjy) 

T33  =  (V  +  ca)(p/2  +  pvrjv/(2c)) 
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T34  =  W  +  ca)/(2crjy)  -  (V  -  ca)/(2c2)  -  rjx2/ ny 
T41  =  V(u2  +  v2)/2 

T42  =  (V  +  ca)(u2  +  v2)pi]x/( 4c)  +  (V  +  ca)pct]  x/(2y  -  2)  -  (V  -  ca)pvrjx/(2rjy)  +  pV 
T43  =  (V  +  ca)(u2  +  v2)p/j  y/(4c)  +  (V  +  ca)pv/2  +  (V  +  ca)pcr{ y/(2y  -  2) 

=  -(V  -  ca)(u2  +  v2)/(4c2)  -  (V  +  ca)v/(2cf]  )  +  (V  +  ca)/(2y  -  2)  +  ur)x  -  vr?x2/r?y 

The  corresponding  equations  for  a  boundary  along  a  ^constant  gridline  have  also 
been  derived  in  this  work.  Those  results  will  be  published  in  the  future. 

The  solution  of  the  boundary  equations  was  made  using  the  same  algorithm  as  the 
interior  equations.  The  numerical  algorithm  is  a  block  tridiagonal  implicit  method  using 
Beam-Warming  second-order  time  differencing. 

RESULTS 

Several  test  cases  were  run  to  test  this  formulation.  The  first  test  case  was  a  solution 
over  a  Cartesian  grid,  where  the  curvilinear  coordinate  formulation  reduces  to  the  boundary 
conditions  given  by  Thompson.1  A  simple  one-dimensional  pressure  pulse  was  given  as 
the  initial  condition,  with  a  freestream  Mach  number  equal  to  0.1  (this  insured  that  there 
was  outflow  only  at  the  downstream  boundary).  A  solution  using  extrapolation  boundary 
conditions  was  obtained  over  a  large  domain,  and  the  solution  was  stopped  before  the 
pressure  pulse  reached  the  downstream  boundary.  Thus,  there  should  not  have  been  any 
influence  from  that  boundary  condition  in  the  interior.  Then,  a  solution  using  the 
nonreflecting  boundary  conditions  was  obtained  on  a  shortened  domain,  and  the  solution 
was  run  until  the  pressure  wave  passed  through  the  boundary.  A  comparison  of  the  two 
solutions  gave  an  indication  of  how  "transparent"  the  nonreflecting  boundary  conditions 
were. 

The  results  of  this  test  were  quite  satisfactory.  There  was  little  noticeable  difference 
between  the  solution  over  the  truncated  domain  and  the  full  solution  (Figure  1).  The  test 
was  repeated  with  a  two-dimensional  pressure  pulse.  In  this  case,  the  results  were  not  as 
good.  Although  the  pressure  wave  did  pass  through  the  boundary,  there  was  some 
distortion  near  the  boundary  (Figure  2).  It  is  possible  that  if  the  cross-terms  in  the 
characteristic  equations  were  not  neglected  in  the  formulation,  the  differences  between  the 
one-dimensional  and  two-dimensional  problems  would  be  insignificant. 

These  tests  were  repeated  on  a  curved  grid  to  test  the  curvilinear  formulation.  These 
tests  gave  the  same  results  as  those  on  the  Cartesian  grid  (Figure  3). 

The  next  test  case  was  that  of  a  single  vortex  passing  through  a  boundary.  The  initial 
condition  was  that  of  a  vortex  velocity  distribution,  but  with  a  uniform  pressure  and  density 
field.  The  low  pressure  within  the  vortex  developed  shortly  after  the  solution  was  begun. 
Again,  both  a  solution  over  a  long  domain  was  found,  with  the  calculation  stopping  before 
the  vortex  hit  the  downstream  boundary,  and  a  solution  using  nonreflecting  boundary 
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conditions  was  found  over  a  shortened  domain.  Solutions  over  both  a  Cartesian  grid  and  a 
curvilinear  grid  gave  the  same  results.  The  solution  over  the  curvilinear  grid  is  shown  here. 
The  freestream  has  a  Mach  number  of  0.2,  and  contours  of  pressure,  vorticity,  and  Mach 
number  are  plotted  (Figures  4-6).  These  results  show  that  the  vorticity  and  the  Mach 
number  propagate  through  the  boundary  with  reasonable  accuracy.  However,  the  pressure 
field  does  not  pass  through  the  boundary  accurately.  Instead,  a  distinct  pressure  reflection 
is  seen,  which  grows  with  time.  Since  the  nonreflecting  boundary  conditions  removes  the 
upstream  influence  of  the  part  of  the  vortex  which  is  downstream  of  the  boundary,  it  could 
be  expected  that  the  pressure  at  the  boundary  would  be  inaccurate. 

The  final  test  case  was  subsonic  vortex  shedding  behind  a  circular  cylinder.  The 
freestream  was  at  Mach  =  0.2,  and  the  Reynolds  number  was  100.  In  this  case,  a  steady 
stream  of  vortices  passes  through  the  downstream  boundary.  The  lift  coefficient  history  is 
oscillatory,  as  shown  in  Figure  7.  At  the  beginning  of  the  calculation,  the  lift  on  the  cylinder 
was  zero,  since  a  symmetric  separation  bubble  behind  the  cylinder  was  formed.  It  was 
necessary  to  perturb  the  solution  to  initiate  the  asymmetric  vortex  shedding.  Both  the 
global  and  the  truncated  solutions  were  run  independently  of  one  another  but  were  started 
in  identical  fashions. 

Once  the  oscillatory  lift  behavior  was  established,  the  two  solutions  differed  in  both 
frequency  and  amplitude.  The  nonreflecting  boundary  conditions,  by  removing  the 
upstream  effects  of  the  vortex  stream  once  it  passed  through  the  boundary,  allowed  a  lower 
frequency  oscillation  than  the  full  solution.  The  amplitude  of  the  lift  coefficient  in  the 
nonreflecting  case  also  seemed  to  level  off  before  the  full  solution  reached  a  constant 
amplitude.  In  the  nonreflecting  case,  a  high  pressure  wave  was  generated  at  the 
downstream  boundary,  which  increased  the  wavelength  of  the  vortex  stream.  Stream 
function  contours  of  the  two  solutions  are  shown  in  Figure  8.  The  elongated  wake  structure 
is  evident  in  the  truncated  solution. 

CONCLUSIONS 

At  present,  nonreflecting  boundary  conditions  present  one  of  the  few  options  for 
treating  outflow  boundary  conditions.  However,  as  the  theory  is  valid  only  for  hyperbolic 
equations,  it  is  necessary  to  understand  the  inaccuracies  that  would  result  when  applying 
them  to  subsonic  viscous  flows.  In  this  work,  it  was  shown  that  errors  in  the  pressure  at  the 
boundary  are  generated,  and  these  errors  propagate  into  the  interior  solution  with  time.  In 
some  cases,  such  as  vortex  shedding  behind  a  circular  cylinder  at  low  Reynolds  number,  the 
outflow  boundary  conditions  eventually  affect  the  entire  wake  structure  and  the  pressure 
forces  acting  on  the  body. 
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Two-dimensional  pressure  pulse  solution.  Initial  condition  P/P 
—  1.  +  0.4  exp(-(x-10)2  -(y-20)2).  a)  nonreflecting  boundary 
condition  applied  at  x=20,  b)  global  solution. 


Figure  3.-  Two-dimensional  pressure  pulse  solution,  same  as  in  Figure  2  but 
on  a  curved  grid.  a)  nonreflecting  boundary  condition,  b) 
global  solution. 
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Figure  7.-  Lift  coefficient  history  for  circular  cylinder.  a)  nonreflectmg 
boundary  conditions  applied  6  diameters  downstream,  b)  global 
solution  with  boundary  16  diameters  downstream.  (Re  =  100, 
M=0.2) 
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APPENDIX  C 


Outflow  Boundary  Conditions  for 
Vortex  Shedding  Behind  a  Circular  Cylinder 


The  flow  about  a  circular  cylinder  was  computed  in  order  to  study  the  effect  of  outflow 
boundary  conditions  on  the  vortex  shedding  process.  Since  the  wake  for  this  flow  persists 
for  a  large  distance  downstream,  it  is  of  interest  to  see  if  a  boundary  condition  applied  to  a 
truncated  domain  can  accurately  reproduce  the  global  flowfield.  If  such  a  boundary 
condition  can  be  found,  the  computational  efficiency  of  these  calculations  will  increase 
greatly,  as  much  smaller  flowfield  domains  may  be  used. 

The  first  step  was  to  compute  the  global  flowfield  to  use  as  a  comparison  to  the 
truncated  solution.  The  numerical  procedure  is  described  here.  A  two-dimensional  finite 
difference  Navier-Stokes  code,  ARC2D,  was  used,  and  the  full  viscous  terms  were  included. 
The  method  used  was  an  implicit  second-order  time-accurate  Beam-Warming  algorithm. 
Second-order  central  differencing  of  the  spatial  derivatives  was  used  with  the  approximate 
factorization  method.  The  calculation  was  laminar,  with  a  Reynolds  number  of  100  and  a 
Mach  number  of  0  1  The  grid  for  the  global  solution  extended  twenty  cylinder  diameters 
downstream. 

The  boundary  conditions  for  the  global  calculation  were  the  default  airfoil  boundary 
conditions  incorporated  into  the  ARC2D  flow  code.1  These  boundary  conditions  are 
described  here. 

For  the  farfield  inflow,  as  well  as  the  top  and  bottom  boundaries  of  the  flow  domain, 
the  following  formulation  based  on  the  one-dimensional  Riemann  invariants  was  used. 
Extrapolation  was  used  to  determine  whether  the  boundary  was  at  an  inflow  or  an  outflow. 
For  an  inflow,  three  variables  were  specified  (Riemann  invariant  Rl,  the  tangential  velocity, 
and  the  entropy)  while  one  variable  was  extrapolated  from  the  interior  (Riemann  invariant 
R2).  For  an  outflow,  R2  was  specified,  and  the  quantities  Rl,  the  tangential  velocity,  and  the 
entropy  were  extrapolated. 

The  Riemann  invariants  Rl  and  R2  are  given  by: 

Rl  =  Qn  -  2.*a/(y-l.) 


and 


R2  =  Qn  +  2*a/(y-U 
with  the  entropy  given  by  p  V/p, 

where  Qn  is  the  velocity  normal  to  the  boundary,  a  is  the  speed  of  sound,  p  is  the  density, 
and  p  is  the  pressure.  For  this  problem,  the  specified  velocities  were  u  =  0.1  and  v  =  0.0 
(based  on  speed  of  sound).  The  freestream  pressure  is  given  by  1  /  y,  and  the  speed  of  sound 
is  1.0. 


On  the  cylinder  surface,  the  velocities  were  set  equal  to  zero,  and  the  pressure  is  found 
from  the  normal  momentum  equation.  An  adiabatic  wall  was  assumed,  and  the  density  was 
extrapolated  from  the  interior.  For  this  global  calculation,  the  outflow  boundary  conditions 
were  all  found  from  an  extrapolation  from  the  interior.  These  boundary  conditions  were 
applied  at  x  =  20. 

The  grid  used  for  this  calculation  was  a  C-type  mesh,  with  dimensions  245x40.  The 
grid  extended  from  -5  <  x  <  20  and  -5  <  y  <  5.  The  cylinder  had  unit  diameter  and  was 
centered  at  the  origin.  The  minimum  normal  spacing  for  the  grid  was  0.016. 

The  truncated  solution  used  the  identical  grid  up  to  x  =  4.  For  this  solution, 
nonreflecting  boundary  conditions  were  used.2  As  the  grid  used  in  this  work  was 
curvilinear,  the  analysis  for  Cartesian  coordinates  given  in  Reference  2  was  rederived  in  this 
work  for  curvilinear  grids.  The  following  results  are  included  here:  both  the  full  and 
truncated  grids,  the  time  history  of  lift  coefficient  for  both  the  global  and  truncated 
solutions,  and  contour  plots  comparing  the  full  and  truncated  solutions. 

Some  comments  on  the  results  should  be  made.  It  is  important  to  note  that  both  the 
global  and  truncated  solutions  were  computed  independently  from  one  another;  that  is,  the 
truncated  solution  was  not  restarted  from  the  full  solution.  Both  solutions  were  started 
semi-impulsively  from  zero.  The  solution  assumes  a  uniform  stream  as  an  initial  condition, 
and  the  cylinder  suddenly  appears  at  time  =  0.  However,  the  code  allows  a  slow  start-up 
time  to  facilitate  the  calculation,  whereby  the  velocity  on  the  cylinder  surface  is  gradually 
reduced  from  freestream  levels  to  zero  over  the  first  twelve  iterations.  Thus,  the  solution  is 
only  "semi-impulsive."  In  addition,  the  code  was  run  with  a  steady-state  version  of  the 
algorithm  to  get  the  solution  started,  and  then  it  was  switched  to  time-accurate.  It  was 
necessary  to  perturb  the  flow  in  order  to  initiate  the  vortex  shedding  process  (seen  clearly  in 
the  CL  vs.  time  curves).  Large  glitches  in  the  solution  may  be  seen  at  y  =  0,  x  =  4  in  the 
truncated  results.  These  are  due  to  a  programming  error  resulting  from  attempts  to 
implicitly  deal  with  an  explicit  wake  cut.  These  errors  have  been  fixed  since  the  time  of 
these  plots. 

The  following  figures  were  taken  from  slides  presented  at  the  Minisymposium  on 
Outflow  Boundary  Conditions  held  in  association  with  the  6th  International  Conference  on 
Numerical  Methods  in  Laminar  and  Turbulent  Flows,  Swansea,  Wales,  U.K.,  on  July  10, 
1989. 
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APPENDIX  D 


Outflow  Boundary  Conditions  Using 
Duhamel's  Equation 


INTRODUCTION 

A  problem  that  has  existed  in  computational  fluid  dynamics  (CFD)  since  its  inception  is 
in  selecting  the  correct  boundary  conditions  for  an  outflow  boundary.  An  outflow 
boundary  is  one  through  which  the  flow  exits  the  computational  domain.  The  errors 
introduced  by  an  inadequate  boundary  condition  can  range  from  a  progressively  inaccurate 
but  physically  realistic  solution  to  a  nonphysical  solution.  The  nonphysical  solution  can 
arise  because  the  boundary  condition  represents  the  effect  of  the  unknown  on  flow  outside 
the  computational  domain,  and  if  the  "exterior"  flow  is  nonphysical,  then  the  "interior"  flow 
is  nonphysical.  A  simple  example  of  such  a  situation  occurs  when  a  steady  state  asymptotic 
boundary  condition  is  used  for  an  unsteady  problem.  In  this  case  all  unsteadiness  in  the 
interior  flow  will  vanish  before  the  boundary  or  the  frequency  of  the  unsteadiness  will 
adjust  such  that  the  boundary  is  a  node  in  wave  space;  either  of  these  situations  is 
nonphysical.  In  general,  the  errors  introduced  by  the  outflow  boundary  condition  are 
greatest  for  unsteady  flows. 

The  early  attempts  at  finding  suitable  outflow  boundary  conditions  concentrated  on 
using  asymptotic  solutions  of  simpler  flow  models,  which  could  be  obtained  analytically. 
For  example,  analytic  solutions  of  the  Prandtl-Glauert  equation  were  used  as  boundary 
conditions  for  the  transonic  small  disturbance  equation.  However,  the  usefulness  of  this 
idea  is  limited  because  of  the  difficulty  in  finding  analytic  solutions  that  will  model  the 
outflow  for  the  Navier-Stokes  or  Euler  equations. 

An  alternative  to  asymptotic  solutions  is  the  idea  of  nonreflecting  boundary  conditions. 
The  physical  basis  of  the  idea  is  that  outgoing  waves  cannot  be  reflected  at  an  (artificial) 
computational  boundary.  A  one-dimensional  analysis  was  developed  by  Enquist  and 
Majda1  for  the  unsteady  Euler  equations,  and  a  version  for  the  unsteady  transonic  small 
disturbance  equation  was  developed  by  Kwak.2  The  ideas  have  been  extended  to  two 
dimensions  by  Thompson,3  but  the  theory  is  really  a  composite  of  two  one-dimensional 
solutions.  Rodman4  has  extended  this  idea  to  curvilinear  coordinates,  but  again  the  two- 
dimensionality  is  "decoupled." 

If  a  particular  flow  is  started  from  zero,  then  the  flow  in  the  exterior  domain  is  exact 
until  waves  from  the  disturbance  in  the  interior  domain  cross  the  boundary.  If  true 
nonreflecting  boundary  conditions  could  be  devised,  then  the  outflow  boundary  condition 
would  be  correct  if  the  computation  was  started  at  zero  time  and  is  time  accurate.  An  error 
in  the  boundary  conditions  due  to  approximations  in  the  formulation  can  give  rise  to  a 
nonphysical  flow  in  the  exterior  domain  and  hence,  ultimately,  the  interior  domain.  A 
principal  difficulty  is  in  finding  an  analytic  relationship  that  will  give  the  boundary 
conditions  for  the  Navier-Stokes  equations.  The  progress  in  nonreflecting  boundary 
conditions  shows  that  finding  a  true  two-dimensional  boundary  condition  is  beyond  the 
present  state  of  the  art. 
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The  domain  to  the  right  of  AB  is  the  exterior  domain,  that  to  the  left  of  CD,  the  interior 
domain.  There  is  an  overlap  of  width  one  grid  cell  in  the  x  direction. 

The  object  is  to  find  the  boundary  condition  on  CD  based  on  information  on  AB.  Note 
that  the  flow  in  the  exterior  domain  is  defined  solely  by  the  correct  data  on  AB. 

If  the  flow  in  the  exterior  domain  can  be  represented  by  a  time  linearized  equation  with 
a  boundary  condition,  U-(t)  on  AB,  where  U.  is  constant  over  the  jth  grid  cell  in  the  y 
direction,  then  the  value  or  Uk(t)  on  CD  at  the  kf)h  grid  cell  is  given  by  Duhamel's  equation 
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where  Ujj^(t)  is  the  "indicial  response"  of  U  on  CD  at  yk ,  that  is,  the  transient  due  to  a  unit 
step  change  in  the  boundary  conditions  on  AB  at  y ■,  lb(t).  £j(t)  is  the  actual  value  of  Lb(t) 
from  the  interior  section.  The  indicial  response  is  obtained  by  running  the  CFD  code  to  be 
used  in  the  analysis  over  a  large  domain  with  a  step  change  in  Un  at  time  zero  and  some  y 
station  and 

U;(t)  =  0  i  *  j  (2) 

where  Lb  (t)  denotes  a  value  at  ys .  The  time  transient  on  the  grid  points  on  CD,  yk  ,  are 
stored  and,  when  divided  by  the  magnitude  of  the  step,  becomes  the  indicial  response 
U|  ^  (t).  In  principle,  this  should  be  computed  for  all  the  j  stations  on  AB  and  constitutes  a 
formidable  mass  of  information.  However,  certain  simplifications  can  be  made.  First,  since 
there  is  no  origin  on  the  infinite  boundary,  AB,  all  that  matters  is  the  location  of  the  yk  point 
relative  to  the  location  of  the  step  change,  that  is,  the  difference  yk  -  y ■,  and  hence  the 
indicial  response  for  the  jth  point  on  A B  is  the  same  as  for  any  other  cell  provided  the 
difference  yk  -  y^  and  the  cell  width  is  held  constant.  The  cell  used  to  impose  the  step 
change  may  be  thought  of  as  the  base  cell.  This  removes  the  amount  of  stored  data  by  a 
factor  of  N,  the  number  of  y  grid  points,  if  a  variable  grid  is  used  then  the  effects  of  varying 
grid  size  on  Ujk  )  can  be  constructed  by  a  sum  of  the  indicial  responses  of  the  base  cell. 
Secondly,  the  time  variation  can  be  approximated  by  curve  fits  which  work  quite  well.  5 
This  removes  the  necessity  to  store  every  time  step;  only  the  coefficients  of  the 
approximating  series  need  be  stored.  In  the  examples  shown  here  this  second  simplification 
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In  Figure  1  the  lift  coefficient  variation  for  the  first  cycle  for  the  transonic  case, 
Mot  =  0.8,  is  shown;  the  boundary  is  0.25  chords  downstream  of  the  airfoil.  The  standard 
code  for  the  truncated  grid  diverges  at  150  time  steps.  The  reason  is  indicated  in  Figure  2 
where  the  pressure  distributions  are  shown.  It  can  be  seen  that  the  nonreflecting  boundary 
conditions  are  totally  inadequate. 

CONCLUDING  REMARKS 


A  new  treatment  of  outflow  boundary  conditions  is  outlined  and  some  preliminary 
results  are  shown.  The  technique  is  applicable  to  a  wide  range  of  equations,  including  the 
Navier-Stokes  equations. 
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APPENDIX  E 


A  Characterization  and  Search  Technique  for 
Unsteady  Flow  Control  Problems 


One  objective  of  the  proposed  work  was  to  use  an  optimizer  to  determine  the  best 
boundary  conditions  for  controlling  a  flow  field.  Mathematical  optimizers  were  available 
for  this  purpose,  but  they  have  several  disadvantages  which  made  them  ill-suited  to  this 
problem.  Consequently,  an  alternative  method  was  developed  which  could  characterize  the 
solution  that  was  being  optimized,  so  that  information  about  all  aspects  of  the  solution 
would  be  found.  For  instance,  various  solution  trends  could  be  documented  as  they  were 
discovered,  providing  information  that  would  help  in  the  search  for  a  new  optimum  if  the 
problem  criteria  were  changed. 

Mathematical  optimizers,  such  as  CONMIN,1  require  the  definition  of  an  objective 
function  which  is  to  be  minimized  by  the  code.  Then,  the  inputs  upon  which  the  objective 
function  depends  are  varied,  starting  from  an  initial  position,  and  an  optimum  combination 
of  the  inputs  is  found.  One  usual  drawback  to  these  optimizer  codes  is  that  only  the  closest 
local  minimum  of  the  objective  function  will  be  found,  therefore  requiring  several  restarts  of 
the  code  with  different  initial  inputs  if  the  problem  is  suspected  to  have  more  than  one  local 
minimum.  Another  disadvantage  is  that  all  quantities  of  interest  must  be  represented  by 
one  function  to  be  optimized,  resulting  in  the  loss  of  information  about  how  they  are 
affected  individually  by  the  inputs.  A  third  drawback  to  these  optimizers  is  that  the  answer 
alone  is  given,  without  significant  other  information,  such  as  how  each  individual  input 
affected  the  answer,  and  if  there  were  any  regions  of  unusual  or  rapid  change  in  the  solution 
field. 


The  current  technique  was  developed  to  address  some  of  these  disadvantages.  It  is 
used  as  a  postprocessor  to  CFD  results.  A  number  of  inputs  to  the  problem  may  be  used, 
and  rather  than  one  objective  function,  several  quantities  of  interest  may  be  specified 
independently.  A  representative  sample  of  inputs  is  used  with  the  CFD  code,  and  the 
corresponding  quantities  of  interest  are  generated  as  outputs.  The  postprocessor  for  this 
information  has  programmed  into  it  the  heuristic  rules  that  a  person  would  follow  if  solving 
the  same  problem.  The  change  in  the  flow  solution  as  a  function  of  one  of  the  variables, 
with  the  other  variables  remaining  fixed,  is  found  for  all  combinations  of  the  variables.  The 
effect  of  one  variable  on  the  solution  may  be  found  to  be  independent,  mildly  dependent,  or 
strongly  dependent  upon  the  other  variables.  In  addition,  how  fast  the  solution  changes 
with  a  given  variable  is  found.  These  computed  slopes  indicate  whether  a  dependence  is 
linear  and/or  monotonic.  The  relationships  between  the  flow  solution  and  the  input 
variables,  along  with  any  regions  of  high  nonlinearity  or  non-monotonicity,  are  documented 
by  the  code.  These  results  give  insight  into  the  cause-and-effect  relationships  between  the 
various  flow  variables  and  computed  flow  field. 

The  postprocessor  may  be  run  iteratively  with  the  CFD  solver  to  find  an  optimum  for 
each  flow  quantity,  or  it  may  be  used  as  a  characterization  code  alone  if  an  optimization  is 
not  required.  For  optimization,  two  techniques  are  used  to  guide  the  search.  First,  a  basic 
hill-climbing  technique  is  used,  where  the  region  surrounding  the  current  optimum  is 
investigated  for  incremental  improvements.  Additionally,  any  exceptional  areas,  as 


shown  in  the  plots  to  indicate  where  this  happens.  In  Figure  2,  a  double  bar  is  shown 
between  the  first  and  fifth  plot,  where  the  function  changes  from  "other"  to  "convex." 

The  code  then  determines  where  future  data  points  should  be  taken,  both  to  search  for 
an  optimum  and  to  investigate  potentially  interesting  regions  with  more  detail.  The 
suggested  values  for  future  iterations  are  based  on  three  criteria.  First,  the  program  looks  in 
the  vicinity  of  the  current  optimum.  It  uses  information  about  the  slopes  and  the  clrms 
value  in  that  area  to  indicate  how  big  a  jump  the  next  iteration  should  take.  The  code  checks 
to  see  if  the  slope  in  that  area  is  much  higher  than  the  median  slope  throughout  the  solution 
field.  It  also  checks  to  see  if  the  optimum  value  found  is  much  higher  than  the  median  value 
throughout  the  solution  field.  If  the  local  slope  or  clrms  value  is  not  significantly  larger  than 
the  median  slope  or  clrms  value,  then  a  relatively  large  increment  will  be  chosen  for  the  next 
iteration.  If  the  local  slope  and  clrms  value  are  high,  then  a  smaller  increment  will  be 
chosen,  especially  if  the  code  determines  that  the  next  iteration  will  fall  in  uncharted 
territory.  If  the  next  iteration  would  be  in  an  area  previously  investigated,  no  new  move  is 
indicated. 

The  second  criteria  for  future  iteration  is  based  on  unusual  slopes  in  the  solution  field. 
If  a  particular  slope  is  much  higher  than  the  adjacent  slopes,  or  is  much  higher  than  the 
median  slope,  a  new  point  is  suggested  that  is  in  the  direction  of  most  improvement.  The 
third  criteria  is  based  on  mismatched  functions  between  adjacent  regions.  If  they  exist, 
points  between  the  two  regions  are  suggested  for  future  iteration.  In  Figures  2-4, 
suggestions  for  future  iterations  are  plotted  on  the  horizontal  axes.  A  square  indicates  a 
search  about  the  local  optimum,  a  circle  indicates  a  search  based  on  unusual  slopes,  and  a 
cross  indicates  a  search  based  on  mismatched  functions. 

Two  iterations  were  performed  after  the  initial  twenty  runs,  and  the  results  are  shown 
in  Figures  5-7.  With  this  much  data,  it  can  be  seen  that  there  are  regions  where  clrms  is  very 
sensitive  to  the  input  variables  and  that  more  resolution  is  necessary  to  understand  the 
results.  The  code  is  programmed  to  flag  any  regions  where  it  is  necessary  to  increase  the 
resolution.  If  a  strong  dependence  on  either  q2  or  q3  within  one  region  shown  in  Figure  5  is 
found,  then  that  region  should  be  subdivided,  so  that  each  resulting  subdivision  is 
independent  of  q2  and  q3.  In  this  example,  the  code  found  that  the  region  where  ql  is  high, 
q2  is  low,  and  q3  is  low  needs  more  resolution.  That  one  region  is  then  subdivided,  and  the 
results  are  shown  in  Figures  8-10. 

The  most  favorable  combination  of  inputs  for  reducing  the  rms  lift  coefficient  was  a 
low  suction  magnitude  with  short  duration  and  interval  pulses.  Figure  11  shows  the  lift 
coefficient  history  for  this  case.  Compared  with  the  no-slip  boundary  condition  case  shown 
in  Figure  1,  the  rms  lift  coefficient  was  decreased  by  10%,  from  0.178  to  0.159,  and  the 
average  lift  coefficient  was  increased  by  5%,  from  0.977  to  1.022. 

This  method  also  investigated  the  region  where  the  solution  was  most  sensitive  to  the 
inputs,  even  though  that  region  did  not  correspond  to  the  optimum  values  found.  The 
characterization  showed  that  the  rms  lift  coefficient  is  most  sensitive  to  long  blowing 
durations,  short  intervals  between  blowing  pulses,  and  high  suction  values.  In  that  region, 
the  most  interesting  result  was  found  where  the  combinations  of  inputs  gave  a  maximum 
rms  lift  coefficient.  The  lift  coefficient  history  is  shown  in  Figure  12.  Although  the  objective 
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Figure  1.-  Lift  coefficient  history  for  no-slip  boundary  conditions. 
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Figure  3.-  RMS  lift  coefficient  value  vs.  pulse  interval  (q2).  Blowing 
duration  (ql)  and  magnitude  (q3)  are  approximately  constant  for 
each  plot.  Full  domain,  first  iteration. 
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RMS  lift  coefficient  value  vs.  blowing  duration  (ql).  Pulse 
interval  (q2)  and  magnitude  (q3)  are  approximately  constant  for 
each  plot.  Full  domain,  third  iteration. 
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Figure  7.-  RMS  lift  coefficient  value  vs.  blowing  magnitude  (q3).  Blowing 
duration  (ql)  and  pulse  interval  fq2)  are  approximately  constant 
for  each  plot.  Full  domain,  third  iteration. 
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Figure  9.-  RMS  lift  coefficient  value  vs.  pulse  interval  (q2).  Blowing 
duration  (ql)  and  magnitude  (q3)  are  approximately  constant  for 
each  plot.  Partial  domain,  third  iteration. 


Figure  11.-  Lift  coefficient  history  for  pulsed  suction  boundary  conditions. 

Suction  =  -0.05*M_. 
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Figure  12.-  Lift  coefficient  history  for  pulsed  suction  boundary  conditions 
Suction  =  -0.174*M ... 


